{
 "metadata": {},
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rank and nullspace of a matrix\n",
      "==============================\n",
      "\n",
      "The following module, rank\\_nullspace.py, provides the functions rank()\n",
      "and nullspace(). (Note that !NumPy already provides the function\n",
      "matrix\\_rank(); the function given here allows an absolute tolerance to\n",
      "be specified along with a relative tolerance.)\n",
      "\n",
      "**rank\\_nullspace.py**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#!python\n",
      "import numpy as np\n",
      "from numpy.linalg import svd\n",
      "\n",
      "\n",
      "def rank(A, atol=1e-13, rtol=0):\n",
      "    \"\"\"Estimate the rank (i.e. the dimension of the nullspace) of a matrix.\n",
      "\n",
      "    The algorithm used by this function is based on the singular value\n",
      "    decomposition of `A`.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    A : ndarray\n",
      "        A should be at most 2-D.  A 1-D array with length n will be treated\n",
      "        as a 2-D with shape (1, n)\n",
      "    atol : float\n",
      "        The absolute tolerance for a zero singular value.  Singular values\n",
      "        smaller than `atol` are considered to be zero.\n",
      "    rtol : float\n",
      "        The relative tolerance.  Singular values less than rtol*smax are\n",
      "        considered to be zero, where smax is the largest singular value.\n",
      "\n",
      "    If both `atol` and `rtol` are positive, the combined tolerance is the\n",
      "    maximum of the two; that is::\n",
      "        tol = max(atol, rtol * smax)\n",
      "    Singular values smaller than `tol` are considered to be zero.\n",
      "\n",
      "    Return value\n",
      "    ------------\n",
      "    r : int\n",
      "        The estimated rank of the matrix.\n",
      "\n",
      "    See also\n",
      "    --------\n",
      "    numpy.linalg.matrix_rank\n",
      "        matrix_rank is basically the same as this function, but it does not\n",
      "        provide the option of the absolute tolerance.\n",
      "    \"\"\"\n",
      "\n",
      "    A = np.atleast_2d(A)\n",
      "    s = svd(A, compute_uv=False)\n",
      "    tol = max(atol, rtol * s[0])\n",
      "    rank = int((s >= tol).sum())\n",
      "    return rank\n",
      "\n",
      "\n",
      "def nullspace(A, atol=1e-13, rtol=0):\n",
      "    \"\"\"Compute an approximate basis for the nullspace of A.\n",
      "\n",
      "    The algorithm used by this function is based on the singular value\n",
      "    decomposition of `A`.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    A : ndarray\n",
      "        A should be at most 2-D.  A 1-D array with length k will be treated\n",
      "        as a 2-D with shape (1, k)\n",
      "    atol : float\n",
      "        The absolute tolerance for a zero singular value.  Singular values\n",
      "        smaller than `atol` are considered to be zero.\n",
      "    rtol : float\n",
      "        The relative tolerance.  Singular values less than rtol*smax are\n",
      "        considered to be zero, where smax is the largest singular value.\n",
      "\n",
      "    If both `atol` and `rtol` are positive, the combined tolerance is the\n",
      "    maximum of the two; that is::\n",
      "        tol = max(atol, rtol * smax)\n",
      "    Singular values smaller than `tol` are considered to be zero.\n",
      "\n",
      "    Return value\n",
      "    ------------\n",
      "    ns : ndarray\n",
      "        If `A` is an array with shape (m, k), then `ns` will be an array\n",
      "        with shape (k, n), where n is the estimated dimension of the\n",
      "        nullspace of `A`.  The columns of `ns` are a basis for the\n",
      "        nullspace; each element in numpy.dot(A, ns) will be approximately\n",
      "        zero.\n",
      "    \"\"\"\n",
      "\n",
      "    A = np.atleast_2d(A)\n",
      "    u, s, vh = svd(A)\n",
      "    tol = max(atol, rtol * s[0])\n",
      "    nnz = (s >= tol).sum()\n",
      "    ns = vh[nnz:].conj().T\n",
      "    return ns"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here's a demonstration script."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#!python\n",
      "import numpy as np\n",
      "\n",
      "from rank_nullspace import rank, nullspace\n",
      "\n",
      "\n",
      "def checkit(a):\n",
      "    print \"a:\"\n",
      "    print a\n",
      "    r = rank(a)\n",
      "    print \"rank is\", r\n",
      "    ns = nullspace(a)\n",
      "    print \"nullspace:\"\n",
      "    print ns\n",
      "    if ns.size > 0:\n",
      "        res = np.abs(np.dot(a, ns)).max()\n",
      "        print \"max residual is\", res    \n",
      "\n",
      "\n",
      "print \"-\"*25\n",
      "\n",
      "a = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])\n",
      "checkit(a)\n",
      "\n",
      "print \"-\"*25\n",
      "\n",
      "a = np.array([[0.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])\n",
      "checkit(a)\n",
      "\n",
      "print \"-\"*25\n",
      "\n",
      "a = np.array([[0.0, 1.0, 2.0, 4.0], [1.0, 2.0, 3.0, 4.0]])\n",
      "checkit(a)\n",
      "\n",
      "print \"-\"*25\n",
      "\n",
      "a = np.array([[1.0,   1.0j,   2.0+2.0j],\n",
      "              [1.0j, -1.0,   -2.0+2.0j],\n",
      "              [0.5,   0.5j,   1.0+1.0j]])\n",
      "checkit(a)\n",
      "\n",
      "print \"-\"*25"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And here is the output of the script."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "-------------------------\n",
      "a:\n",
      "[[ 1.  2.  3.]\n",
      " [ 4.  5.  6.]\n",
      " [ 7.  8.  9.]]\n",
      "rank is 2\n",
      "nullspace:\n",
      "[[-0.40824829]\n",
      " [ 0.81649658]\n",
      " [-0.40824829]]\n",
      "max residual is 4.4408920985e-16\n",
      "-------------------------\n",
      "a:\n",
      "[[ 0.  2.  3.]\n",
      " [ 4.  5.  6.]\n",
      " [ 7.  8.  9.]]\n",
      "rank is 3\n",
      "nullspace:\n",
      "[]\n",
      "-------------------------\n",
      "a:\n",
      "[[ 0.  1.  2.  4.]\n",
      " [ 1.  2.  3.  4.]]\n",
      "rank is 2\n",
      "nullspace:\n",
      "[[-0.48666474 -0.61177492]\n",
      " [-0.27946883  0.76717915]\n",
      " [ 0.76613356 -0.15540423]\n",
      " [-0.31319957 -0.11409267]]\n",
      "max residual is 3.88578058619e-16\n",
      "-------------------------\n",
      "a:\n",
      "[[ 1.0+0.j   0.0+1.j   2.0+2.j ]\n",
      " [ 0.0+1.j  -1.0+0.j  -2.0+2.j ]\n",
      " [ 0.5+0.j   0.0+0.5j  1.0+1.j ]]\n",
      "rank is 1\n",
      "nullspace:\n",
      "[[ 0.00000000-0.j         -0.94868330-0.j        ]\n",
      " [ 0.13333333+0.93333333j  0.00000000-0.10540926j]\n",
      " [ 0.20000000-0.26666667j  0.21081851-0.21081851j]]\n",
      "max residual is 1.04295984227e-15\n",
      "-------------------------"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}